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Abstract 

We present a theoretical model for Bragg scattering from a Bose-Einstein conden- 
sate (BEC) in the vicinity of a magnetic Feshbach resonance, using a two c-field for- 
malism, one c-field for the atom and the other for a molecule formed of two atoms. 
We use this model to numerically simulate a recent experiment |Tj investigating the 
effects of strong interactions on the Bragg spectrum from a ^^Rb BEC. Results from 
these simulations are in very good quantitative agreement with the experimental re- 
sults, confirming the importance of the resonance bound state in the dynamics of the 
condensate for fast experiments like Bragg scattering. 

1 Introduction 

Bose-Einstein condensates (BEC) v^^ith weak interparticle interactions have been, in many 
cases, successfully described using a pseudopotential formulation, in which just two pa- 
rameters are involved: the density and the s-wave scattering length. Furthermore, in a 
large proportion of situations, mean-field theory can be used, leading to a description in 
terms of the Gross-Pitaevskii equation. Arising from this success, a quest for a tunable 
and possibly large interaction strength began, leading to the study of systems in which the 
scattering length was the result of a Feshbach resonance, whose use made it possible to 
tune the interatomic interaction of a BEC over a wide range. 

Using Feshbach resonances, it has become possible to investigate condensates with 
strong interparticle interactions. This has been done both experimentally and theoreti- 
cally with quite wide success, even though there are reasons to question the validity of the 
standard theoretical procedures at some of the interaction strengths used. The pseudopo- 
tential and mean-field theory methods are the result of a perturbation treatment, which 
must definitely fail for sufficiently high interaction strengths and densities. 

With this in mind, in Paper I (2] and Paper II [3] we introduced a more careful treat- 
ment of interactions mediated by weakly bound molecular states, such as arise in a Fesh- 
bach resonance. This was done by introducing a molecular field, whose interaction con- 
stants can be determined phenomenologically from scattering length and binding energy 
data. In our treatment, the interaction constants are relatively weak, but nevertheless 
reproduce many of the results of a simple pseudopotential method, especially for static 
properties, such as the condensate shape. Building on this, in Paper II we formulated a 



Bogoliubov description, and showed that there are changes in the excitation spectrum at 
higher energies. Using this description, we found modifications in the Bragg scattering 
spectrum from a homogeneous infinite condensate very similar to those which were ex- 
perimentally found by Papp etal. |T] 

What this means is that, by treating the dynamics behind the change in interaction 
strength in the Feshbach resonance, we find that a mean-field treatment of a strongly in- 
teracting system is still very much applicable, provided that there are two mean fields, one 
for the atoms and one for the molecules. 

In this paper we continue the work presented in our two previous papers, and inves- 
tigate a realistic system. We apply the formalism of Paper I and Paper II to the specific 
case of a recent experiment by Papp etal. [T], in which the excitation spectrum of a Bose- 
Einstein condensate of ^^Rb was measured using Bragg scattering, near the Feshbach res- 
onance at 155 G. This experiment was deliberately designed to explore a region of parame- 
ter space in which the perturbation theory would not be expected to be valid. And indeed, 
by tuning the scattering length to large values, they found significant deviations from the 
Bragg scattering behaviour predicted by the simplest perturbative and mean-field theo- 
ries. 

The results of our work are very satisfactory. Although the experiment was not de- 
signed to test this kind of theory, and thus some significant parameters are hard to es- 
timate, we obtain quantitative agreement with their experimental results, with no fitted 
parameters. 

1 . 1 Properties of the Bragg spectrum 

The excitation spectrum of a homogeneous condensate for large momentum transfer is 
given by the sum of the kinetic energy and the chemical potential of the condensate. 



where k is the photon momentum, n is the density of the condensate and as is the s-wave 
scattering length. In the case of an inhomogeneous condensate, for example a condensate 
in a harmonic trapping potential, it was found by Stenger et al. [4] that this formula can be 
used provided n is interpreted as the density-weighted density of the condensate. In their 
work, the use of the density-weighted density was theoretically justified by using a local 
density description of the condensate. Blakie et al. |5 6| simulated the system using the 
Gross-Pitaevskii equation, and confirmed the basic validity of this approximation. In the 
region where (TJ is valid, it is equivalent to the Bogoliubov excitation spectrum in the limit 
of large A;. 

The prediction (1) is expected to be valid as long as the condensate is dilute [na^ « 1), 
the excitation is in the free-particle regime m » 1, where ^ - [Snnas)'^'^ is the conden- 
sate healing length) and the scattering amplitude is momentum independent {kag « 1). 
The aim of the experiment of Ref fTl was to investigate the properties of the Bragg spec- 
trum in a region where the scattering length is large. This means that the condensate in- 



teractions cannot be treated as mean-field (y 8nna^ ~ 0.5), the excitations are not clearly 
particle-like (A;^ ~ 2), and the scattering amplitude is not clearly momentum independent 
(fcas~0.8). 

1.2 Experimental results and issues 

In the experiment of Ref. (J, the shifts of the Bragg spectra for large scattering lengths 
showed a significant deviation from the theoretical predictions based on {1). The ex- 
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perimental results were also compared to theoretical predictions outlined in detail in jT], 
which are not in agreement with the experimental data. 

Another theoretical model is presented by Kinnunen et a/. |8] , who studied Bragg spec- 
troscopy from a uniform, strongly interacting ^^Rb condensate using time-dependent Hartree- 
Fock-Bogoliubov theory. They took into account the momentum dependent scattering 
amplitude, but found only qualitative agreement with the experiment. 

1 .3 Interpreting the Experiment 

There are several issues that complicate the analysis of the experiment of Ref. [Ij . We will 
address the most important of these in the following: 

1. Initial state : The nonlinear effects that were investigated in the experiment are 
more pronounced the less dilute the condensate is. The density of the initial state in 
the experiment is therefore enhanced by a series of ramps of the scattering length. 
Creating an initial state in this way makes the experimental procedure even more 
complex and takes the system further away from the ideal case studied theoretically 
by □ and [8). 

Furthermore, the initial state parameters are not explicitly defined, which makes 
analysis of the line shift result difficult, since this depends on the properties of the 
condensate at the onset of the Bragg pulse. 

2. Inhomogeneity : The trapped condensate is spatially inhomogeneous, and also strongly 
time varying, because of both three-body losses and condensate expansion in the 
breathing modes, and even as a result of the Bragg scattering process itself (as we 
shall see in section|4T). Papp etal. measure and use space averaged densities, rather 
than the density weighted densities, which (as we noted above) are more appropri- 
ate when comparing with results for a homogeneous condensate. In addition to this, 
they also average densities over the duration of the experiment. 

The problem with using the space-averaged density is that unless the condensate 
has a clearly defined volume, the space-averaged density cannot be accurately de- 
termined. In an experiment such as that of 1 1 1, the volume is not easily determined 
and has to be approximated in one way or another. 

In [T], the time- and space-averaged density was determined by assuming that the 
density profile of the condensate is given by a Thomas-Fermi profile with a width 
given by a variational solution to the Gross-Pitaevskii equation (GPE). The varia- 
tional model, outlined in |9| is, however, not an accurate representation of an exact 
solution of the GPE, which even in three dimensions, is not very difficult to find nu- 
merically. 

Furthermore, it is not clear how accurate a description of the shape of the conden- 
sate is given in this case by a Thomas-Fermi profile. As we shall see in section fOI in 
our simulations the shape of the condensate is very different from that given by the 
Thomas-Fermi approximation, both before and during the application of the Bragg 
pulse. 

3. Variable pulse length and intensity : The condensate density varies more rapidly 
in time as a result of Bragg scattering at larger scattering lengths. By introducing 
the condition that the density of the condensate cannot change by more than 30% 
during the Bragg pulse, based on predictions from the variational model, the exper- 
iment is forced to use progressively shorter Bragg pulses for larger scattering length. 
To make sure that roughly the same quantity is scattered out each time, the intensity 
of the pulses is appropriately increased. 
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Scattering length flj [ool 
Figure 1: Frequency scales in the experiment [l]. 



In our calculations we find that the processes involved are sensitive to the duration 
and intensity of the Bragg pulse, because the condensate expands, because there 
are three-body losses, and because the spectroscopic resolution improves for longer 
pulses. It is therefore important to reproduce the experimental parameters as faith- 
fully as possible. Unfortunately, neither duration nor the intensity of the Bragg pulse 
are explicitly stated in [IJ, so we have inferred their values from the spectra and the 
limits on the number of Bragg-scattered atoms. 

4. Time scales : Bragg scattering is a fast process and it is therefore important to con- 
sider the other time scales associated with the experiment; if other processes occur 
on a time scale similar to that of the Bragg scattering, it is likely that those processes 
are important for the dynamics of the condensate. 

There are three time scales that are relevant in this type of experiment, shown in 
Fig-ffl in terms of their corresponding frequencies: The frequency of the applied 
Bragg pulse, the inverse of the pulse duration and the binding frequency of the 
bound state in the Feshbach resonance. For large scattering lengths and short Bragg 
pulses, these frequencies are comparable, and it is therefore very probable that the 
bound state dynamics become important to the overall dynamics of the experiment. 



2 Formalism 

The formalism of Paper I and Paper II proceeds in brief as follows: To model the Feshbach 
resonance bound state, we add an additional field, corresponding to a bound atom pair 
(refered to as a "molecule"), to the usual Hamiltonian for a trapped system of interacting 
Bosons. The equations of motion for the atom field y/ and molecule field cp in the resulting 
c-field model are given by 

ifi-^ = -——yfix] + Ta{Vaix]yf{x)+Uaamx)fv^{x) + gy/*{x)cl){x)} 
at 2m 

-ij[\^fr{x)f + 2\if>{x)ffy/{x), (2) 



i?i^^ = -- — (l){x) + ymUe+V,nix])cf>{x) + ^y/\x]\, (3) 
at 4m ^ 2 > 
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where Uaa - 4nfi?ai,gl m is the background interaction strength, Vq and Vm are the ex- 
ternal trapping potential for the atoms and molecules respectively. The last term in j2) is 
added to account for losses from the condensate due to three-body recombination events 
[TOl ; we discuss this more extensively in Sect. 13.2.3] 



2.1 Projectors 

Ta and T m are the atom and molecule projectors that restrict the wavefunctions to the low 
energy subspace below the momentum cutoff, 
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where is the Heaviside step function, and the momentum space cutoff in the j-th direc- 
tion is given by 



max(A;j) 



(6) 



These projectors arise from the simulation grid, and they are defined in order to avoid 
aliasing in our simulations [TT|. However, as discussed in Paper I, it is necessary to include 
the momentum space cutoff for two other reasons: in order for a pseudopotential treat- 
ment of the interaction to be valid and in order to avoid excessive quantities of the noise 
being added in the c-field method. 

2.2 Parameters 

The parameters g and e are the coupling strength and detuning in the Feshbach resonance 
respectively. As shown in Paper I, they are in our formalism given by 



e = 



f^2^2 (;r-2Aa,)(l-^f(f)) 



2m 



Aasa+tif))-n 



2Aas[\+t[f,]-n 



(8) 



where t{x) - x - arctanl/jc, and h^a^/m is the molecular binding energy corresponding 
to the s-wave scattering length as I2j . The parameter A is the renormaUzation factor, de- 
termined by the momentum space cutoffs. 

The effect of the Bragg field on the condensate is included by making the following 
substitutions in the equations of motion: 
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where 



Vopt - Vqcos[x- q - (t)t), 



(10) 



where q and cj are the wavevector and the frequency of the Bragg pulse respectively f5lF6j , 
and Vq is the amplitude of the optical potential, given in terms of the Rabi frequency O 
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and the excited state detuning A, 




(11) 



The optical potential for the molecule is chosen to be twice that of the atom on the as- 
sumption that the atoms in the molecule are very weakly bound, and for these purposes 
behave almost independently. 

2.3 Renormalization factor 

We tend to refer to the parameter A as the momentum space cutoff. However, in reality 
this is only true in the special case of isotropic cutoffs. The relationship between the renor- 
malization constant A and the momentum space cutoffs kx,cut, fcy,cut and fc^^cut. is given 



where V is the ellipsoidal volume spanned by the momentum space vectors, correspond- 
ing to the projectors (4][5). 

In the simplest case, the momentum space cutoff is the same in all directions, and the 
volume of the populated low energy subspace is spherical so that evaluating (12) gives 



where fc^^cut is the value of the isotropic cutoff. In the case of an anisotropic cutoff, as 
in the case for the numerical calculations in this paper, the exact value of A needs to be 
evaluated using equation i fTZt . More details are given in Appendix lB.2l 

3 Simulations 

We simulate the experiment of Ref. fT\ by numerically solving the equations of motion 
(2) and l|3) in three dimensions. To model the effects of quantum fluctuations in c-field 
theory, the wavefunctions in (2) and (3) have a random amplitude added to the initial 
states, corresponding to half a virtual particle per mode [11] . 

3.1 Momentum space truncation 

The trapping potential in the experiment of [U is cigar-shaped, with an aspect ratio of 
1 : 46. This, along with the fact that the Bragg pulse is applied in the axial direction, and the 
Bragg momentum is relatively large, leads to a system that is computationally demanding. 
To include all the relevant physics, and at the same time ensuring that the c-field methods 
are still valid, and that the system is still computationally tractable, we make a truncation 
of the momentum space, neglecting all the modes that do not make a significant contribu- 
tion to the dynamics of the system. This procedure, which involves dividing momentum 
space into bands, each centered around one of the Bragg orders, expresses the wavefunc- 
tions y/ and cp as 
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where y/„ and (pn are the Fourier transforms of the momentum space wavefunction in the 
band centred around riQ. The details of the procedure are given in Appendbc lB.ll 

We find that only the four momentum bands corresponding to the orders n- -1,0,1,2 
acquire sufficient population to affect the simulation. For each band there wiU be an atom 
wavefunction i//„ and a molecule wavefunction . The equations of motion for the atom 
wavefunctions are given by 

at 2m [ 2 

+ Uaa [Aof-i + A-ilj/Q + A-zfi + A-3f2) + g (Vo0-1 + Vl^Po + V20l) 

-ijlCoyz-i + C-iy/Q + C-2V1 + C-3f2]} (16) 

= -^—^Vo + J'a\Vay/o + — 
at 2m { 2 

-iY[Ciyj-i + CoVo + C-11//1 + C-2V2)} (17) 

or 2m { 2 ^ I 

+ Uaa [A2f-\ + Alfa + Aq^i + A-1^2] 

+g{vl^2 + Vl<Pi + V*-i<Po) - ir{Coy/i + CiVo + C2V-1)} (18) 
= -^^:^V2 + ya\vaV2 + ^Vie~''^' + Uaa[A3yr_i + A2y/o + Aiy/i + Aoy/2] 



n + Val VaVfo + Y (v-1 + Vie"-"] 



dt 2m 

+g{V-i(f>i + Vo(t>2) - ir{C3f-i + C2y/o + C11//1 + Coy/2)} , (19) 
where for brevity we have suppressed the spatial and temporal dependence, and where 

Vl^V^ + i2nQ—-n^Q^, (20) 
dx 

and the factors A,,, B,i and C„ are given inAppendbc lB.il 

Similarly, the equations of motion for the molecule wavefimctions become 

if,^!tl ^ -^-^(^_i + T,Jy,„(^_i + Vo(^oe'"" + gi/'-iVo} (21) 
of 4m ^ ' 

= -^cPo + 'J'rn\v„,<l)o + VoU-ie-'''' + cp,e""] 
at 4m ' ^ 

+ |(2i/._ii//i+V'o)} (22) 



at 4m I V ; 

+g{y/-iy/2 + VoVi)} (23) 

,^d^ ^ -^02 + 3'm{V;,7(/)2 + ^^o(/)ie-''"* + f (2v'oV2 + V'D}- 
at 4m I 2 ' 



(24) 



The projectors 7a and Vm are defined as in i|4l[5) with the fc^;- directional cutoff now 
given by 

Afc 

fcx,cut=— , (25) 
4 

where AA; is the width of each momentum space band. The projectors are the same for 
all the wavefunctions iffn and (/>„; for each band the projectors are ellipsoids centered 
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Table 1: Experimental average densities (column 3) for different scattering lengths (column 1). The 
data for the shift (column 2) are taken from Fig. 3(a) in [IJ. Column 4 lists the approximate ranges 
of the average densities, where the values in brackets are our estimates. Column 5 shows the time 
averages of the density- weighted density that we use in our numerical calculations. 



as 


/shift 


Time- and space-average 


Range 


Time-average density-weighted 


[flol 


[kHz] 


density [10^^ m"^] 


[10^^ m"^] 


density [10^^ m"^] 


150 


0.9 


7.6 


7.6 


10.857 


300 


1.7 


7.2 


7.2 


10.286 


500 


2.9 


7.3 


7.1-7.6 


10.429 


585 


3.0 


6.5 


(5.6-7.4) 


9.286 


695 


3.8 


6.9 


(6.2 - 7.6) 


9.857 


805 


3.9 


6.1 


4.9-7.4 


8.714 


890 


4.6 


6.5 


(4.8 - 8.2) 


9.286 



around the midpoint ky - kz - Q and kx - nQ. The full wavefunctions, rfr and (p, are 
thus projected onto four disjoint regions in momentum space (see Fig.lSlin Appendix IB. It . 
The band width Afc is chosen as a compromise between two factors: it needs to be large 
enough to include as much as possible of the momentum space wavefunction, but at the 
same time small enough to not include too much of the initial noise. It is also important 
that the individual bands are not overlapping. 

3.2 Simulation parameters 

In choosing the parameters for our simulation we follow the experimental setup as closely 
as possible, and model a condensate of 40,000 ^^Rb atoms in a trap with cylindrical sym- 
metry and an aspect ratio of 46.2 [Vz - 2.9Hz, Vr - 134Hz). The scattering length as ranges 
from 150ao to 890^^o, and we use relationship between the scattering length and the pa- 
rameters g and £ derived in Paper I and given by equations (7) and l|8) . 

3.2.1 Initial state 

In the experiment an initial condensate was created with a scattering length of ISOaq • The 
scattering length was then ramped to an unspecified low value, exciting the large ampli- 
tude breathing modes in the condensate. At the inner radial turning point of the breathing 
mode oscillation, the scattering length was ramped up to the desired value, and the Bragg 
pulse was applied. Through this process the condensate becomes much denser, making 
the nonlinear effects on the Bragg spectra more clearly visible. However, it is hard to know 
exactly what the initial state at the commencement of the Bragg pulse is; had the conden- 
sate not been compressed in this way, the initial state would have been clearly defined. 

In our simulations we create the initial state for the Bragg spectroscopy by performing 
the following steps: 

1. We set the scattering length to a small, arbitrary value flinit, typically of the order of 
a few flo- 

2. We numerically solve the time-independent equivalents of the equations of motion 
(2) and (3) for this scattering length, using the Thomas-Fermi equations that we de- 
rived in Paper II as the starting position. 

3. We quickly ramp the scattering length up to the value of interest, in the range be- 
tween 150(20 and SSOao- The speed of the ramp never exceeds dslus < 0.25hl ma^. 
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Table 2: Properties of the Bragg pulse used in our simulations for different scattering lengtlis. The 
values of the duration are inferred from the spectral widths in ilj . 



Scattering length as 

[flol 



Duration t 

[ms] 



Amplitude Vq 
[h kHz] 



Scattered fraction 



150 
300 
500 
585 
695 
805 
890 



0.09 



0.45 
0.24 
0.14 



0.11 
0.10 



0.12 



0.80 
0.94 



0.13 
0.40 
0.67 



0.80 



1.07 



6.2% 
6.7% 
6.5% 
7.5% 
6.4% 
7.1% 
8.0% 



4. We then apply the Bragg pulse while continuing to run the simulation, and calculate 
the resulting time-average of the density weighted density for the duration of the 
Bragg pulse. In our simulations, the density weighted density h{t) is given by 



Here N{t) and h{t) involve "effective numbers of atoms", counting each molecule as 
two atoms, corresonding to what would in practice be measured in an experiment. 

5. We adjust flinit appropriately, and redo steps 1-4 until the time-averaged density 
weighted density obtained matches that of the experiment. 

The experimental space- and time-averaged density is inferred from the predictions 
of the line shift (Fig. 3(a) in yj). This varies from 7.6 x lO^^cm"^ for 150ao to between 
approximately 4.8 and 8.2 x lO^^cm"^ for 890ao, see Tab.[l] We relate the experimental 
space-averaged density to our density-weighted density by noting that the space-averaged 
density for a Thomas-Fermi profile is 0.7 times the density averaged density. The same 
factor is not necessarily right for other profile shapes, and as we shall see, the condensate 
in our simulations is quite far from Thomas-Fermi shaped. However, we believe that this 
nonetheless gives us the best estimate of the density used in |l] that we can reasonably 
expect to get, since it corresponds to the procedure used in the experiment to estimate the 
space-averaged density. 

3.2.2 Bragg pulse 

The Bragg pulse, modelled by fTO) and assumed to be square, is applied at the start of the 
simulation with a wavenumber of A; = 47r/780 nm in the axial direction of the condensate. 
The pulse durations are not explicitly stated in [T|, but can be inferred from the data for 
the widths of the spectra (Fig. 3(b) in pT]), where the contribution from the pulse duration 
will be inversely proportional to the pulse length as Aw - 0.36/ ?puise (where the number 
0.36 comes from the rms width of a Gaussian fit to the Fourier transform of a square Bragg 
pulse), see Tab.|4]in AppendixjA] 

Thus, the pulse duration, and therefore also the simulation time, ranges from 0.09 ms 
for 890(20 to 0.45 ms for ISOaq. see Tab.[2l In Tab.|2] we have also listed the Bragg pulse 
amplitudes Vq for the different values of the scattering length. The intensities of the Bragg 
pulse are not stated in fll, but as in the experiment we have chosen it so that we always 
have between 5% and 10% of the condensate being scattered, see Tab.m 




(27) 



(26) 
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3.2.3 Three-body loss 



The rate of particle loss from the condensate arising from three-body recombination events 
varies with the scattering length, approximately proportionally to a* p2), and is also ex- 
tremely sensitive to the density. There are no exact values of the three-body loss coefficient 
(7 in 12)) available; here v^^e have used the theoretical values given by Braaten et al. (13), 
which qualitatively agreed with previous experimental data from Roberts et al. [14) . How- 
ever, as is noted in 1 13| these values are highly uncertain. 

The way in which to include three-body loss in a c-field formalism was originally de- 
veloped by Norrie et al. (10] . This treatment includes a stochastic term in the equations 
of motion, but Norrie shows that this term can in most cases be neglected to a good or- 
der of approximation. In (2), we have used this approximate form, and have extended it 
phenomenologically to include the molecule population as well. This extension is quite 
simple minded. We include all of the losses in the equation for the atomic field, and use the 
total density of atoms plus molecules in the loss term as the object that corresponds most 
logically to the measurable density of atoms. Since the atomic field is very much larger 
than the molecular field in the situations we are considering, this kind of model should be 
a reasonably accurate approximation. 

In principle, the formulation of the theory in terms of atoms and molecules provides an 
opportunity to give a model of three-body loss which would incorporate the actual mech- 
anism of three-body loss as arising from inelastic collisions between atoms and molecules. 
In such a collision, both the atom and the molecule would normally be transferred to un- 
trapped states, and be lost from the system — equivalent to a loss of three atoms. We hope 
to develop this kind of model in a future publication. 



3.3 Simulations of structureless atoms 

For comparison with the simulations based on our formalism, we also run simulations 
based on a simple GPE. This corresponds to modeling the condensate using a single field 
^ , and by letting the interaction strength be determined solely by the scattering length. In 
the c-field formalism, the equation of motion for the single-component condensate is in 
this case given by 

= >I'W + T„{VaW>I'W + Lfol^(A:)|2vp(x)}-i7|¥(x)|4¥(x), (28) 

at 2m 

where the parameters are the same as in Equation (2), except the atom-atom interaction 
which is now given by 

t/o=— — (29) 
mil-2Aas/7i) 



4 Results of Simulations of the Mean- Field Equations 

The underlying equations of motion in the c-field formalism are the same as those of mean 
field theory, and quantisation is introduced by the inclusion of fluctuations in the initial 
state. The inclusion of the fluctuations can cause very dramatic changes in the nature 
of the solutions, as was found in [15] . In the Bragg scattering problem under study here, 
we have found that the effects of the quantum fluctuations are in fact rather small. It is 
therefore logical to study first the solutions of the equations in the absence of the added 
noise in the initial conditions, which amounts to a mean-field description of the system 
of atoms and molecules. Indeed we find that these simulations provide a very satisfactory 
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Figure 2: Coordinate space profiles for flj = 890flo at a time (a) just before the Bragg pulse is applied, 
and (b) at the end of the Bragg pulse. The left panels show a slice of the condensate in the axial 
direction, i.e. the direction of the pulse, and the right panels show slices in the radial direction. 
In the bottom right panel we show two slices, corresponding to the crest (solid line) and trough 
(dashed line) of the interference fringe centered around x = 0. The inset in the bottom left panel 
shows the profile appearance in the region marked by the red box. The parameter xq is the length 
scale associated with the x-axis of the trap, given by xq = \/hl2m(j)x ~ 6.665 x 10~^m. 



description of the problem, which agrees very well with the experimental results of [T). 
The effect of the noise terms is thus a matter of determining relatively small corrections to 
the mean field theory, and this will be done in the following section. 



4.1 General behaviour 

The coordinate space profiles from a typical simulation run are shown in Fig.|2] The scat- 
tering length is in this case 890flo- Fig. |2(a)| shows the radial and axial profiles of the con- 
densate after it has been ramped to the scattering length of interest, at the moment just 
before the Bragg pulse is applied. Fig. |2(b)| shows the the profiles for the same simulation 
run, at the end of the Bragg pulse. 

As can be clearly seen in Fig. |2(a]| the condensate profile in the axial direction before 
the onset of the pulse is similar in shape to a Thomas-Fermi profile, whereas in the radial 
direction it is more like a Gaussian. This is a result of the elongated shape of the conden- 
sate, due to the aspect ratio of the trap. As noted in Sect ll.3l this makes any estimate of the 
average density difficult to justify. 

The Bragg pulse is applied in the axial direction; in Fig. |2(b]| we can clearly see the effect 
of this as interference fringes in the axial profile of the condensate. In the radial direction 
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Table 3: Condensate three-body loss and density change during the Bragg pulse, for the different 
scattering lengths. 



Scattering length as [flol 


Three-body loss 


Density change 


150 


1% 


70% 


300 


2.5% 


50% 


500 


2.5% 


33% 


585 


4% 


30% 


695 


5% 


30% 


805 


5.5% 


26% 


890 


6% 


24% 



we have therefore plotted two distinctly different profiles, corresponding to the crest and 
trough of the central fringe. We obtain similar profiles with large density variations for 
each of the different scattering lengths in our simulations. 

The particle losses arising from three-body recombination events and the change in 
density in our simulations are very different from those predicted by fTl , whose predic- 
tion is that the density will change by "less than 30%". In contrast, in our simulations the 
density changes by up to 70% of the initial density (see Tab.[3l) Furthermore, in the ex- 
periment, the three-body loss is observed to be "typically <30%" [I]; whereas in our sim- 
ulations the losses never exceed 10% of the total atom number. We believe that the main 
reason for these differences is the inappropriate model used to describe the condensate 
in 111. However, there is also a significant degree of uncertainty in our calculations of the 
three-body loss, because of the lack of accurate data for the loss rate, and this could also 
contribute to the discrepancy. 



4.2 Bragg spectra and lineshift 



The Bragg spectrum is obtained by changing the frequency difference o) in fTol , and calcu- 
lating the momentum transferred to the condensate for each frequency. We calculate the 
normalized momentum transfer as 



P{t) = j^^^^ f dk{\y/{k, r)|2 + 2\ip{k, f)|2) k. 



(30) 




15 20 
Bragg frequency [kHz] 



Figure 3: The Bragg spectra for three different values of the scattering length: 150ao (blue triangles), 
500flo (red circles) and 890flo (black squares). The solid lines are Gaussian fits to the data points. 
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— Experimental data 
-■ — Simple shift prediction 

Atom Bogoliubov 
H — Atom-molecule Bogoliubov 

— * Atom simulation 

-0 Atom-molecule simulation 
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Scattering length as [ao] 
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Figure 4: The shift of the peak of the Bragg spectra for different values of the scattering length. The 
solid lines are results from our simulations, the dashed lines are the experimental results and predic- 
tions presented in Q] and the dash-dotted lines are calculations based on the Bogoliubov treatment 
of Paper II. Our atom-molecule simulation (green diamonds) is significantly different from that of 
the structureless atom (red stars), but agrees well with both the experimental data (black circles) 
and our atom-molecule Bogoliubov calculation (magenta plusses). The prediction of the lineshift 
based on the excitation spectrum in the large k limit of equation (T] (blue squares) shows a very dif- 
ferent behaviour for large scattering lengths, as does the structureless Bogoliubov calculation (cyan 
crosses). The error bars on our atom-molecule calculation indicate the uncertainty in the exper- 
imental estimates of the density, as is shown in Tab.[T] similar error estimates would apply to the 
other curves, but have been omitted for clarity. 



Typical spectra from our simulations are shown in Fig.|3l The difference in width be- 
tween the different spectra is due to the change in duration of the Bragg pulse. The peaks 
of the Bragg spectra in Fig.|3]are shifted from the position of the corresponding peaks for 
the non-interacting gas, located at approximately 15.4 kHz. 

In Fig.lUthe shifts of the Bragg spectra from result for the noninteracting case are plot- 
ted as a function of the scattering length. 

For comparison, we have also plotted the experimental results of the Bragg lineshift 
from Ref. \lj . Fig.|4]also includes the theoretical prediction of the lineshift based on (1) , but 
using the density weighted density instead of the space-averaged density used in the ex- 
perimental paper. For a Thomas-Fermi profile the density weighted density is 10/7 times 
larger than the space-averaged density, as discussed in Sect. l3.2.n This correction elimi- 
nates the anomaly apparent in Fig. 3(a) in Q], in which the experimental data and the sim- 
ple shift prediction agree almost perfectly up to a scattering length of about 500ao. and 
then deviate sharply. Using the density- averaged density, a smooth increase in deviation 
is apparent. 

As can be seen clearly in the figure, the simulations based on our model show quan- 
titative agreement with the experimental data. We have included error bars on the ex- 
perimental data points; these indicate the uncertainty in the experimental estimates of 
the density in the experiment. Similar error bars should therefore also be included on the 
other lines in Fig.lH but we have omitted these for clarity. 

For small scattering lengths, the molecule field is very small, but neverthless plays an 
important role since its presence gives rise to a positive scattering length, in contrast to 
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the negative background scattering length. The fact that the binding energy is larger at 
low scattering lengths makes it possible for the molecule field to adiabatically follow the 
atom field, and thus the condensate behaviour is very similar to that predicted by a GPE 
description. This is clear in Fig.lD where we have included the result from the simulations 
based on the GPE (28). At larger scattering lengths, the bound state evolves more slowly 
and the atom-molecule simulations become very different from those for structureless 
atoms. 

Finally, in Fig.|4l we have included the results from Bogoliubov treatments of the ideal 
case of a uniform condensate, both for the case of a single atom field, as in l5j|6J, and 
in the case of an atom-molecule system, as in Paper 11. We find that the atom-molecule 
Bogoliubov treatment shows surprisingly close agreement with both our simulations and 
with the experimental data. 



r(fc) 


>{ 




V2 




V2 



5 Results of Full C-Field Simulations 

In the c-field methods, the effect of quantum fluctuations is included by adding stochastic 
terms to the initial state, corresponding to on average half a particle per mode. In our 
treatment of Bragg scattering, we will follow the approximate procedure as noted in [TT| of 
adding Gaussian random noise, r{k) and s(7c), with zero mean and unit standard deviation 
to the initial momentum amplitudes for the atoms y/Q and (po, according to 

r(k) , . s(7c) 

yf(k)^yfo{k] + —, (/,(7c) = (/,o(fc)-H— . (31) 
V 2 V 2 

Each simulation run can then be seen as corresponding to a single run of an experiment, 
and the expectation values of observables are obtained by taking the average of several 
different runs. The average of the noise amplitudes is obviously 

2. (32) 

corresponding to half a noise particle per mode. 

Fig.EJshows the phase of the spatial atom and molecule fields for a slice in the xy-plane 
for the same system as in Fig.|2]at the end of the Bragg pulse. Since the initial stochastic 
terms are added to the momentum space wavefunctions, and since the molecule projector 
encompasses a much larger part of momentum space than the atom one, there are many 
more noise particles in the molecule field than in the atom one. Despite this, and the fact 
that the molecule field is much smaller than the atom one, there is still a clearly visible 
phase coherence in the molecule field. 

It is remarkable that the noise evident in the phase of the molecule field has very lit- 
de effect on the results of simulations. The large positive scattering length arises directly 
from the population of the molecule field, and one might have expected its value to be 
significandy affected by the quantum fluctuations as they appear in the c-field model. 



5. 1 Density Weighted Density in Terms of C-Fields 

The correct computation of the density-weighted density involves some care, since it in- 
volves products of four field operators, including both molecule and field operators. The 
details of how this is done are presented in Appendix[Cl whose results are in summary: 

1. The average total particle number in the noise simulations is given by 



Nit): 



/ dx\{na{t)}sy^ + 2{n,„[t)}sy^- — -Ar, 



(33) 
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Figure 5: Phase plots ofaslice ofthe condensate centered around z = for a scattering length of flj = 
890flo. showing the atom part (top panel) and the molecule part (bottom panel) of the system. The 
total number of particles is around 40,000, with approximately 34,000 atoms in the form of atoms, 
and 6,000 atoms as molecules. The total number of noise particles is 30,000 for the atom field and 
500,000 for the molecule field. The parameter xq is the length scale associated with the x-axis of the 
trap, given by xq = ^/h/lmwx = 6.665 x 10~^m. 



where {ha(x, ?)}sym and {hm{x, f)}sym are the symmetrically ordered averages 

{naCx.Dlsym = \fHx,t)y/{x,t)\ , (34) 

' I J sym 

{nm{x,t)}sym = W(x,t]<p[x,t)\ . (35) 
' I I sym 



The parameters and Am corresponds to the noise on the atom and molecule co- 
ordinate space wavefunctions, respectively, given by 

(|r(fc)|2) = Aa, (36) 
(Ismf) = Am- (37) 

Since the molecule field has much more initial noise added to it, is much larger 
than Aa- 

2. The coordinate space density-weighted density is given by 

nii) = ==^Jdx[n^(f)-H4n2^(f)-2A,„{na(f)}sym-2Afl{«m(f)}sym 

+4{/la(f)}sym{"m(?)}sym + AflA,„]), (38) 



where ria [x, t) and ix, t) are given by 


nl{x,t) = {nl(.x,t)}^^-2Aa{na{x,t]}^^+^, (39) 

a2 



,{x,t) = {nUx,t)}^^-2Am{nmix,t)]sym + ^, (40) 



15 



11 




7.5 I ' ' ' ' 1 

0.02 0.04 0.06 0.08 

t [ms] 

Figure 6: Time evolution of the density-weighted density during the application of the Bragg pulse 
to a condensate. The scattering length is flj = 890flo. and the pulse length is 0.09ms. The solid 
line is calculated using 1381 for 30 different runs, with error bars indicating the statistical error from 
these runs. The dashed line shows the density-weighted density from a simulation without initial 
noise. The dotted lines shows the time- average of the density- weighted density for the noise free 
simulation (blue) and for the full simulations (magenta). 



where 

{"a(^-f)}sym = {t^^' (x, DV'' (x, f)}^^ , (41) 
{"i(^'«}sym = {4>^\x,t]<pHx,t)]^^. (42) 

We run 30 simulations v^^ith noise and compare the density-v^reighted density obtained 
from these runs using (38) with that obtained in a single simulation run without any stochas- 
tic terms added to the initial state. The result for one of these comparisons is shown in 
Fig.[6l where we have plotted the evolution of the density- weighted density for a scatter- 
ing length of as - 890^0 • As can be seen clearly in the figure, the initial density- weighted 
density is the same for both the noise-free simulation and the average of the 30 runs with 
noise. However, as the condensate evolves over time, the result from the noise simula- 
tions is slightly lower than that from the noise-free run. The resulting time-average of the 
density-weighted density will therefore be slightly higher if we neglect the initial fluctua- 
tions, although the size of the change is much less than the experimental uncertainty. 



5.2 Bragg Spectra from C-Field Simulations 

For the full c-field simulations, instead of the results in (27] [30), we calculate the momen- 
tum transfer as 

= - ([ dki\y/ik,t)\^ + 2\(l)ik,t)\^-^]k), (43) 

N{t)\q\ \J [ 21/ 



where the factor of | is subtracted to account for the initial noise. Similarly, the total num- 
ber of particles is given by 



N{t)^(^j dki^\yr{k,t)f + 2\cl)(.k,t]\^-^j}. (44) 
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Figure 7: Bragg spectrum for a scattering lengtli of 890flo. sfiowing the result obtained from averag- 
ing over 30 runs of simulations performed with stochastic terms added to the initial modes (magenta 
dots), as well as the result from a single simulation without initial noise (blue stars). The solid lines 
are Gaussian fits to the data points. 



Fig.[7]shows the momentum transfer calculated using equation 1431 , where the average has 
been taken over 30 simulation runs. In comparison, we have also plotted the momentum 
transfer from a single simulation run writhout any initial noise terms. The noise simula- 
tions give a spectrum that is slighdy narrower than the noise-free simulation and with a 
slightly lower amplitude. However, the position of the spectral peak is essentially the same 
for both the simulations with and without initial noise. 

In our simulations, we find that averaging over several different noise simulations in 
this way gives us very results similar to those obtained by running the same simulation 
without including the noise. Although the vacuum fluctuations seem to have some small 
effect on the evolution density, overall, the effect on the condensate dynamics appears to 
be unimportant to the Bragg scattering experiment. We can therefore be confident that 
the simulations of the mean-field equations of the atom-molecule system which we did 
in Sect.|4]— equivalent to omitting the initial quantum fluctuations — provide a reliable de- 
scription of the Bragg scattering experiment. 

6 Conclusions 

The aim of the experiment of Ref (Jl was to to probe the behaviour of a Bose-Einstein 
condensate in the regime where two major simplifications normally made in its theoretical 
description were not valid. These simplifications are made in terms of three dimensionless 
parameters: 



1. Weak Interactions : This requires y Snna^ « 1. It is important to note that this ap- 
proximation is necessary not only for the validity the Gross-Pitaevskii equation, but 
also for the validity of the local quantum field theory, to which the Gross-Pitaevskii 
is an approximation. 

In the experiment the condensate was compressed, and the scattering length in- 
creased by using a Feshbach resonance, in order to ensure the violation of this con- 
dition. 

2. Local Interactions : By this, it is meant that the length scale on which processes of 
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interest take place is much larger than that of the interactions. In the experiment the 
momentum transfer involved in the Bragg scattering was chosen to be sufficiently 
large that the momentum dependence of the scattering amplitude would be impor- 
tant. 

In addition, the parameters of the experiment were chosen so that the the relevant quasi- 
particles, that is, those with momentum corresponding to the Bragg wavenumber, were 
definitely not in the free particle regime. 

In the three papers in this series we have shown how to take account of all of these 
within a tractable formalism. The most significant aspect of both the experiment and the 
theory is the clear demonstration that the large scattering lengths generated by Feshbach 
resonances do not give rise to interactions of the hard-sphere kind, as treated originally 
by Huang and Yang I16II17I . Indeed, it is remarkable that the classical Huang- Yang theory 
works so well for systems with Feshbach resonance enhanced interactions. For this rea- 
son, in Paper II we investigated stationary states, the Thomas-Fermi approximation and 
the Bogoliubov excitation spectrum of our model of coupled atoms and molecules, and 
in fact found that even when ag - 890flo. the corrections were quite modest, though quite 
perceptible. The Bragg scattering experiment is essentially a measurement of the excita- 
tion spectrum, and the frequency changes it presents are a measure of the deviation from 
the spectrum expected of the corresponding hard sphere model. These corrections are in 
fact quite modest; only about 10% of the actual Bogoliubov quasiparticle frequencies. 

6. 1 Relating Theory to Experiment 

The experiment set out to test the limits of conventional theory, and convincingly achieved 
that aim. However the procedure used was not ideal for comparison with our detailed 
model. The most challenging problem is the absence of any measurements of the ini- 
tial state of the condensate immediately before applying the Bragg pulse. The issue is 
further complicated by the procedure used to enhance the density of the condensate, be- 
fore ramping the scattering length to the value used for the Bragg pulse. The result is an 
initial state for the Bragg scattering which, not being a stationary state, cannot be defini- 
tively determined. In order to compare our computations with experiment we have relied 
on the time and space averaged density measurements implicit in their presentation of 
the frequency shifts expected from the Huang- Yang theory. We have converted these to 
the appropriate values of the time-averaged density-weighted density, and using these we 
achieve our results, which are in very good quantitative agreement with the measured re- 
sults. 

We would consider it of importance in any future experiments to present 

1. Either : Measurements of the initial state; 

2. Or : A precise quantitative description of the procedure used to create each initial 
state from the initial condensate, which can be reliably modelled as a stationary 
Bose-Einstein condensate. 

The presentation of results as time averaged quantities should be avoided; these create 
very significant computational difficulties. 

6.2 Further Opportunities 

The methods we have developed can clearly be applied to other problems in which the 
flexible adjustment of the scattering length afforded by Feshbach resonances has been ex- 
ploited, for example the Bose-Nova problem, and the related problem of bright solitons. It 
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is also conceivable that the methodology could be extended to study the physics of Efimov 
states in the presence of a Bose-Einstein condensate. 
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A Simulation parameters 

The values of different parameters used in the simulations are listed in this section. The 
duration of the Bragg pulse is inferred from the data for the widths of the Bragg pulse (Fig. 
3(b) in (D), where the contribution from the pulse duration will be inversely proportional 
to the pulse length as Aw - 0.36/ tpuise- The duration for the different scattering lengths 
are listed in Tab.m where we also list the values of the three-body loss parameter y. The 
loss parameter has been determined by using the corresponding values of the three-body 
recombination rate K3 given by Braaten et al. [13) . 

Table 4: Pulse lengths and three-body loss parameters for different scattering lengths. The data for 
the width are taken from Fig. 3 in Q] , and values of the loss parameter are calculated using values of 
the three-body recombination rate in [l3l . 



Scattering length 


Width from duration 


Pulse length 


Three-body loss 


as [flol 


A/ [kHz] 


t [ms] 


7 (cm~^/s) 


150 


0.8 


0.45 


5 X 10"^^ 


300 


1.5 


0.24 


1.5 X 10-27 


500 


2.5 


0.14 


5 X lO-^'' 


585 


3.0 


0.12 


1.5 X 10-26 


700 


3.4 


0.11 


2 X 10-26 


800 


3.7 


0.10 


2.5 X 10-26 


890 


4.0 


0.09 


3 X 10-26 



B Projectors and momentum space truncation 
B. 1 Momentum space truncation 

To include all the physics that we are interested in, the momentum space needs to include 
at least the first order Bragg momentum, at - \q\. Assuming the Bragg pulse is only 
applied in the x-direction, we can write the optical potential as 

l^opt = Vocos[Qx-Mt) = ^(e'W^-<^« + e-im-o't)^^ ^455 

where Q- \ q\. 

To fulfill this condition as well as the condition that the number of grid points Ntextgrid - 
2" for some positive integer n, the number of grid points in the x-direction is chosen to 
be Nx - 2048. The y- and z-directional grids are chosen to have Ny-Nz- 64, since these 
directions are of less importance, making the total number of grid points 2 x 2048 x 64 x 64. 
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This would not only make the simulations very computationally heavy, but also, since in 
the c-field formalism we will have on average half a quasiparticle of noise per mode in the 
initial state, we would get many more noise particles than condensate particles. Accord- 
ing to the validity condition for the c-field methods [10], this would make our simulations 
invalid. 

To get around this problem but stUl include all the momentum space of relevance, 
we neglect the parts of momentum space where the population will be insignificant, and 
include only those modes that are initially populated or where we can expect to get signifi- 
cant population from scattering. Because the interest here is Bragg scattering with a Bragg 
pulse applied in the positive x-direction, we divide momentum space into bands in this 
direction, each centered around nQ for some n, where Q is the momentum of the pulse. 
We can thus write the wavefunctions y/ and (p as 

V^W^^V/nWe""?', (46) 

n 

(/)W = ^(/)„(x)e''"'3^, (47) 

n 

where i//„ and (p„ are the Fourier transforms of the momentum space wavefunction in the 
band centred around riQ. 

The projectors and CPm are given by 



= e 







■ + + ■ 



J.2 1.2 1.2 

V '^x,cat '^y,cut '^z,ciit I 



h.2 J.2 J.2 



(48) 



(49) 



with the x-directional cutoff now given by 
Afc 

fcx,cut=— , (50) 
4 

where Afc is the width of each momentum space band. The projectors are the same for all 
the wavefunctions iff„ and </)„; for each band the projectors are ellipsoids centered around 
the midpoint at y = z = and x - nQ. 



B. 1 . 1 Four significant bands 

We find that only the four momentum bands corresponding to the orders n - -1,0, 1,2 
win be significant populated during our simulations. We then have 

f(x) = i//_i(x)e"''^'' + y/oix) + (x) e''?'' + i//2(jc)e^''^'' (51) 
(/)(x) = (X) e'"^'' + (po{x) + (pi [x] e"^'' + (p2 {x] e^''^''. (52) 

Fig.[8] shows the projector for the full wavefunctions for this case in the xy-plane, where 
we have also indicated the width of each band A A; and the momentum space cutoffs fc^.cut 
and fcy.cut- 

This gives us the following expression for the squared norm of i// 

\\ff{x)f = A-3{x)e~''^'' + A-2[x)e~"^'' + A-i{x)e~'^'' + 

+Ao[x) + Ai(x)6?''3^ + AzMe^"^"" (53) 
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Figure 8: Atom field projector (top panel) and molecule field projector (bottom panel) in the xy- 
plane for the case of four bands in momentum space being significantly populated. The width of 
each band is AA;. The blue areas indicate the regions of momentum space that the wavefunctions are 
projected into, determined by the parameters fcx,cut and fcy,cut- The parameter xq is the length scale 
associated with the x-axis of the trap, given by xq = i/ hl2m(Ox ~ 6.665 x 10"'^ m. Here Afc = 1.3x 10~^ 



m, kr 



■- 3.2 X 10" m, and k 



y.cut 



:3.5x 10"^ m 



where 

^0 = |l//-i|2 + |l//o|2 + |l/Ai|2 + |i^2l^ (54) 

Ai = vUVo + vlVi + vlVi^ A*_^ (55) 
A2 = ^ff*_yy/i + y/oy/2^ A*_2 (56) 
A3 = y/liy/2^A*_3 (57) 

Similarly for the squared norm of (p we have 

|(/)(x)|^ = B-3{x)e~''^'' + B-2[x)e~"^'' + B-i[x)e~"^'' + 

+Bo[x) + Bi{x)e'Q'' + B2{x]e^'-^'' (58) 

where 

Bo = \'P-xf + m^ + \(pi\^ + \(p2f (59) 
Bi = (p*_^(Po + iPl(Pi + (PliP2^B*_^ (60) 

B2 = (p*_^tpi + yj*Q(p2^B*_2 (61) 

B3 = (/)!i(/)2 = B!3 (62) 

The density squared now becomes 

[\y/[x)f + \(p[x)\^f = C_3(x)e"^''3-" + C-2(x)e"2''?'^+C-i(x)e"''?^' + 

+ Co {X) + Ci [X) e'^"" + C2 ix) e^'^"" + C3 (jc) e''^'' (63) 
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where 



Co = {Ao + 2Bof + 2\Ai+2Bi\^ + 2\A2 + 2B2\^ + 2\A3 + 2B3f (64) 
Ci = 2(Ao + 25o)(Ai+25i) + 2(A_i + 2B_i)(^2 + 2B2) 

+2(^_2 + 2B_2)(^3 + 2B3) = C!i (65) 

C2 = (^i + 2Bi)2 + 2(^o + 2Bo)(^2 + 2B2) + 2(^-i + 2B_i)(^3 + 2B3) = C!2 (66) 

C3 = 2(^0 + 2Bo)(^3 + 2^3) + 2(^1 + 2Bi)(^2 + 2B2) = C!3 (67) 

and terms with n < -3 or n > 3 have been neglected. 

These expressions, together with the expressions for y/ and (p are substituted into the 
equation of motion for the atom wave function l|2). Dropping all terms with a factor of 
ginQx ^j^j^ -1,0, 1,2, and collecting terms corresponding to the same band together, 
we get 

of 2m [ 2 

+ Uaa [Aaf-i + A-ilj/Q + A-zfi + A-3f2) + g [vl(p-i + y/*i(f>o + y/lipi) 
-ij[Cof-i + C-ii/'o + C-2V1 + C-3f2]} (68) 

+ Uaa {Aiy/-i + Aoy/o + A-iy/i + A-zy/z) + g (V-i0-i +fo(po + fl(pi + Vl^Pi] 
-/7(CiV-i + Coi/'o + C-iyji + C-2V2)} (69) 

" -^Vi + 5'a|K,Vi + y(voe-'"^ + i/'2e'"^) 

+ Uaa [A2f-\ + Alfa + Aq^Ji + A-iyj2] 

+g[y/l(p2 + Vl<pi + V-\<Po) - '7 (CoVi + CiVo + C2V-1)} (70) 

ot 2m [2 

+g[V>-i'Pi+V*Q'p2) - ir[C3V-i + C2l/'0 + C11//1 + Cof2)}, (71) 

where 

= V^ + /2«— -n^Q^. (72) 
dx 

Similarly, the equations of motion for the bands of the molecule wavefunction 0, be- 
come 



+ ?a { Van + y {v-ie-""' + Vie""'] 



ih 



= -^^^fl)_^ + y,^\v,n(l)-i + Vo(l)oe''" + gy/-iyfo\ (73) 
of 4m ^ ' 

= -^cl)o + yrn{v,n<l>o + Vo[cP-ie-''"' + cPie''''] 

+ |(2v._iV/i+Vo]} (74) 

b /i^V^ 
ih-^ = ——^cl)i + y,n\Vr„(pi + VoUoe-''"' + cl,2e''"'] 
at 4m I >. 

+g{y/-iy/2 + VoVi)} (75) 

= -^02 + 3'm|v^m(/)2 + l^o(/)ie-'"* + f (2v'oV2 + V'Dl- (76) 
of 4m I 2 ) 
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B.2 Renormalization parameter 



In order to calculate the renormalization constant, we need to evaluate the integral 
r dk 

'-l^^' '''' 
where the relationship between the renormalization factor, A, and the integral is 3 = \nK. 

In the simplest case, the momentum space cutoff is the same in all directions, and 
the volume of the populated low energy subspace is spherical. In this case we can use 
spherical coordinates, to get 

J = ^ = \ndp = 4;rfcB,cut, (78) 

Jo Jo Jo p Jo 

where fcfl,cut is the value of the cutoff, so that we get the simple relationship between the 
renormalization factor and the cutoff A = A;^ cut- 



B.2.1 Anisotropic cutoff 

In the case of an anisotropic momentum space cutoff, calculating the renormalization 
factor becomes slightly more complicated. In our simulations we have, as is described 
in section IrTI a momentum space that is divided into four bands. We therefore have to 
calculate 



n n JVn K 



(79) 



where Vn is the volume of the low energy subspace in band n. 

Each band in the truncation has an ellipsoidal projector, symmetric in the yz-plane, 
with maximum value fcy cut - ^z.cut = fcyz.cut- In the x-direction, each band is centered 
around kx-nQ and they all have the same width Afc. 

Because of the cylindrical symmetry of the volume we can simply the problem by 
changing to polar coordinates to get 



JnQ-&kl2 Jo 

Jo 

where pmax is given by 



nQ+^kl2 />pniax(0 2np 

2 J2^P^^ 



Q-Akl2 Jo P +( 

nQ+Akl2 

2n\ (In (p|,„, «•)+<•')- In (C2))d<: 



Pmax(0 - 



yz, cut 



We then get 



/■Afc/2 / 

„ = 2;rj^ In 



{nQ + Ak/2y 



""yz.cut 



[Akliy 



+ 2nQ(+{nQf + k 



yz.cut 



-21n(C+nQ) 



(80) 
(81) 

(82) 

d(. (83) 



Since the volume ¥„ is an ellipsoid and not a sphere, we have that fc^^^^, ^ (AA:/2)^, 
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and since Afc < 2Q, i.e. the bands are not overlapping, this integral has the solution 



X In 



l/T(A;fi,cut)-l ^'^^ 

'f{2nQ- 2fcfl,cut ^l-T(fcfl,cut)-T(?lQ)] +4(1- T(fcB,cut)] ((?2Q)^ - fc^.cut) 

^f[2nQ + 2fcfl,cut Vl-T(fcW-T(«Q)) +4(1- T(fcB,cut)) [inQ)2 - fc^^^^J 



4n^nO 
+ — In 

1-T(fcii,cut) 

where t(x) =x^/{Ak/2f. 



nQ + Akl2 



- 27inQln [nQ [nQ + Akl2)) , 



(84) 



C Density- weighted density 

We wish to calculate the density- weighted density « ( f ) for our coupled atom and molecule 
system in the Wigner formalism used in this paper. 



C.l Wigner ordering 

For an operator a, we know the symmetrically ordered average, 

^ ( -2 ^f2 -t - -^t ' -t2 -t - ->t 't2 -~2l 

= -<^a a' + a' aa' a+ aa' a+ aa' aa + a a j . 

Assuming that the commutator is 
- A, 

and that 



N- d^d 



we find 



I J svm 2 



) sym 

Since we also have 



we get 



^' = {^f'}sym-2A{N},y„ + y. 

We therefore get the averages 



N 
N2 



sym 2 



(85) 

(86) 
(87) 

(88) 

(89) 

(90) 

(91) 
(92) 
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C.2 Atom-molecule density-weighted density 

We now consider the case of an atom operator f) and a molecule operator (f>{x, t), with 
commutators 

[l/A(x,f),1//'^X,f)] = ^a, (93) 

(p{x,t),(pHx,t]^ = Am, (94) 

The average total atom number for this system is given by 

N{t) = ^ J dx^^Hx,t)ii{x,t) + 2(p\x,t)^{x,t)\^, (95) 

where, as usual, we count a molecule as two atoms. Using the commutation relations this 
can be expressed as 

lV(7) = ^J dx\^{na{x,t)]^^ + 2{n,n{x,t)]^y^-^-^r^, (96) 
where {ha[x, t)]^^^ and {hm{x, f)}syin are the symmetrically ordered averages 

{n«(A;,f)}sym = \vf\x,t)T4f{x,t)\ , (97) 

' I J sym 

{«„,(x,f)}sym = U\x,t)^[X,t)\ . (98) 
I I svm 



The density-weighted density for the system is given by 

•\2 



n{t]^^=(^^ dx\ijj'^{x,t)i'{x,t) + 2(p^(x,t)^{x,t)^ J- (99) 
Using the same approach as in the previous section, we can express this as 

^= (^j dx\^nl ix, t] + 4nl,{x, t) - 2^m{naix, t)]syni - 2Aa{nm {x, r)}sym 



nit) 

mt) 



+4{naix,t)]sy^{nr,jix,t)]sy^ + AaAm)), (100) 



where rig (.x, t) and [x, t) are given by 



[x,t) = {nl{x,t)}^^-2Aa{na(x,t)}^^ + ^, (101) 



n 



A2 



[x,t) = {nf„(x,f)} -2A„,{n,„(x,r)}sym+^, (102) 



where 



{"«(^'«}sym = {v>^'(x,f)V>'(^,fl}^^, (103) 
{"if^'fJlsym = {<P^^ix,t)(p^{x,t)]^ . (104) 



C.3 Check with initial state 

The initial state corresponds to the two states y/{x) and (p{x), given by 
f(x) 

y/ix) = il/o{x) + ——, (105) 
v2 

<^(;c) = cl)o[x) + ^, (106) 
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where 



(|rW|2> 
(|fW|*> 

Then we have 



= A 
= 2A 



2 



sym 



({"flWIsym) = 



) sym 



/ fix) 

A^ 

|i/'oW|* + 2A«|i//oW|^ + -^, 

li/^oWl +— , 
six) 

(l)0ix) + —pr 

V2 



^2 

= |0o(x)|%2A,„|0oW|^ + -^, 



|VoWH(/>oW| 



1 \Voix)\ 



Aak>oW A„A„ 
! !_ + 

2 4 



And therefore 

_^ /£^x[|yoW|^ + 2|(/)oWp]^ 

"~ /dx[|l/AoW|2 + 2|(/)oW|2] ' 

as expected. 



(107) 
(108) 
(109) 
(110) 



(111) 
(112) 

(113) 
(114) 



(115) 



(116) 
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